Etude de marché pour l’exportation de poulet

Objectif de la mission

L’entreprise La poule qui chante souhaite se développer à l’international. L’objectif est d’obtenir une liste de pays cibles pour l’exportation du Poulet.

• Les Indicateurs utilisés : Cette étude de marché va suivre la méthode PESTEL, où pour chacun des pays, un ou plusieurs indicateurs politiques, économiques, sociologiques, technologiques, environnementaux et légaux seront analysés.

  • Politique : Stabilité politique (indice)
  • Economique : PIB par habitant ($); Taux de croissance (PIB par habitant, % annuel); Investissements étrangers directs (% du PIB)
  • Sociologique : Population totale, Population active (%)
  • Technologique : Accès Electricité (%), Performance logistique (indice)
  • Environnemental : Emission de CO2 (t)
  • Légal : Solidité des garanties juridiques (indice)

De plus, le niveau de consommation (t) , d’exportation (t), et d’importation (t) de poulet pour chacun des pays sera pris en compte.

Les jeux de données proviennent de deux sources principales : l’Organisation des Nations Unies pour l’Alimentation et l’Agriculture (FAO) et La Banque Mondiale.

• Postulat de départ : Les pays cibles seront ceux qui importent une quantité non négligeable de poulet et ayant des bons indicateurs socio-économiques.

• Méthode : Afin de trouver les meilleurs pays cible qui répondent aux critères ci-dessus,nous effectuerons une analyse sur les variables par la méthode d ’Analyse en Composantes Principales (ACP) pour réduire le nombre de variable à prendre en compte. Ensuite, nous réaliserons un regroupement des pays par la méthode de la Classification Ascendante Hiérarchique et par la méthode des kmeans sur les axes de l’ACP puis sur les variables de départ. Les pays en commun entre ces quatres modèles déterminera la liste définitive.

I) Préparation de l’environnement R et importation des données

#vidage de la mémoire
rm(list = ls())

I.1) Initialisation des packages

library(ggplot2) # Pour les graphes
library(ggrepel) # Pour les graphes
library(cowplot) # Pour les graphes
library(viridisLite) # Pour les graphes
library(dplyr) # Pour les modifs des datas
library(tidyr) # Pour les modifs des datas
library(stringr) # Pour les modifs des datas
library(FactoMineR) # Pour l' ACP et clustering
library(factoextra) # Pour l' ACP et clustering
library(NbClust) # Pour l' ACP et clustering
library(corrplot)
library(ggpubr)
library(gplots) # pour les heatmaps
library(heatmaply) # pour les heatmaps

library(ggrepel)

library(leaflet) #Pour les cartes
library(rgdal) #Pour les cartes

options(scipen=999)

I.2) Importation et préparation des données

Données issues de la FAO

Afin de préparer le jeu de données final, des modifications préliminaires sont nécessaires :

• Modification du nom de la colonne ‘Area.Code..ISO3.’ par ‘ISO3’

• Filte ede certains pays : suppression du Pays China (qui regroupe mainland, hong kong, taïwan etc…) et remplacement du code ISO3 “41” de China mainland par “CHN”.

• Sélection de l’année 2019, où l’on dispose de plus d’informations

• Suppression des colonnes de descriptions pour chacun des jeux de données

• Les modifications spécifiques sont commentées au niveau des lignes concernées

PIB2010_2019 <- read.csv("~/Documents/OpenClassRooms/Projet9/Data/PIB_2010-2019.csv", encoding="UTF-8")
names(PIB2010_2019)[names(PIB2010_2019) == "Area.Code..ISO3."] <- "ISO3"
PIB2010_2019=subset(PIB2010_2019, ISO3 != 'CHN')
PIB2010_2019$ISO3[PIB2010_2019$ISO3 %in% "41"] <- "CHN"
PIB2010_2019 <- subset(PIB2010_2019, Element.Code == '6119') #Nous gardons uniquement le PIB par habitant. 
PIB2010_2019=PIB2010_2019[,c(3,4,5,6,7,8,10,11,12)]
PIB2019 <- subset(PIB2010_2019,  Year == 2019)


PoliticalStability2010_2019<- read.csv("~/Documents/OpenClassRooms/Projet9/Data/Political_stability2010-2019.csv", encoding="UTF-8")
names(PoliticalStability2010_2019)[names(PoliticalStability2010_2019) == "Area.Code..ISO3."] <- "ISO3"
PoliticalStability2010_2019=subset(PoliticalStability2010_2019, ISO3 != 'CHN')
PoliticalStability2010_2019$ISO3[PoliticalStability2010_2019$ISO3 %in% "41"] <- "CHN"
PoliticalStability2010_2019=PoliticalStability2010_2019[,c(3,4,5,6,7,8,10,11,12)]
PoliticalStability2019 <- subset(PoliticalStability2010_2019,  Year == 2019)



Population2010_2019<- read.csv("~/Documents/OpenClassRooms/Projet9/Data/Population_2010-2019.csv", encoding="UTF-8")
names(Population2010_2019)[names(Population2010_2019) == "Area.Code..ISO3."] <- "ISO3"
Population2010_2019=subset(Population2010_2019, ISO3 != 'CHN')
Population2010_2019$ISO3[Population2010_2019$ISO3 %in% "41"] <- "CHN"
Population2010_2019=Population2010_2019[,c(3,4,5,6,7,8,10,11,12)]
Population2010_2019$Value <- Population2010_2019$Value*1000
Population2019 <- subset(Population2010_2019,  Year == 2019)


FB2015_2019 <- read.csv("~/Documents/OpenClassRooms/Projet9/Data/FB2015_2019.csv", encoding="UTF-8")
names(FB2015_2019)[names(FB2015_2019) == "Area.Code..ISO3."] <- "ISO3"
FB2015_2019=subset(FB2015_2019, ISO3 != 'CHN')
FB2015_2019$ISO3[FB2015_2019$ISO3 %in% "41"] <- "CHN"
FB2015_2019=subset(FB2015_2019, !(ISO3 %in% c(248,228,15,51,62,'ANT','SCG'))) #Ces pays qui ne dispose n'existent plus dans l’interval de temps d'intérêt. 
FB2015_2019=FB2015_2019[,c(3,4,5,6,7,8,10,11,12)]
FB2015_2019$Value=FB2015_2019$Value*1000
FB2019 <- subset(FB2015_2019, Year == 2019)

Données issues de la Banque mondiale

Nous analyserons par la suite uniquement les données de 2019. Pour les données issues de la Banque Mondiale, nous garderons uniquement la colonne du code ISO3 des pays et la valeur des données pour l’année 2019

Acces_electricite <- read.csv("~/Documents/OpenClassRooms/Projet9/Data/Acces_electricite.csv", encoding="UTF-8", sep=";")
Acces_electriciteF = Acces_electricite[c(2,9)]
names(Acces_electriciteF)[names(Acces_electriciteF) == "Country.Code"] <- "ISO3"
names(Acces_electriciteF)[names(Acces_electriciteF) == "X2019"] <- "Acces_electricite"


Emission_CO2 <- read.csv("~/Documents/OpenClassRooms/Projet9/Data/Emission_CO2.csv", encoding="UTF-8", sep=";")
Emission_CO2F = Emission_CO2[c(2,9)]
names(Emission_CO2F)[names(Emission_CO2F) == "Country.Code"] <- "ISO3"
names(Emission_CO2F)[names(Emission_CO2F) == "X2019"] <- "Emission_CO2"


Indice_juridique <- read.csv("~/Documents/OpenClassRooms/Projet9/Data/Indice_juridique.csv", encoding="UTF-8", sep=";")
Indice_juridiqueF = Indice_juridique[c(2,11)]
names(Indice_juridiqueF)[names(Indice_juridiqueF) == "Country.Code"] <- "ISO3"
names(Indice_juridiqueF)[names(Indice_juridiqueF) == "X2019"] <- "Indice_juridique"

Indice_perf_logistique <- read.csv("~/Documents/OpenClassRooms/Projet9/Data/Indice_perf_logistique.csv", encoding="UTF-8", sep=";")
Indice_perf_logistiqueF = Indice_perf_logistique[c(2,11)]
names(Indice_perf_logistiqueF)[names(Indice_perf_logistiqueF) == "Country.Code"] <- "ISO3"
names(Indice_perf_logistiqueF)[names(Indice_perf_logistiqueF) == "X2019"] <- "Perf_logistique"

Investissement_etranger <- read.csv("~/Documents/OpenClassRooms/Projet9/Data/Investissement_etranger.csv", encoding="UTF-8", sep=";")
Investissement_etrangerF = Investissement_etranger[c(2,9)]
names(Investissement_etrangerF)[names(Investissement_etrangerF) == "Country.Code"] <- "ISO3"
names(Investissement_etrangerF)[names(Investissement_etrangerF) == "X2019"] <- "Investissement_etranger"

taux_croissance <- read.csv("~/Documents/OpenClassRooms/Projet9/Data/taux_croissance.csv", encoding="UTF-8", sep=";")
taux_croissanceF = taux_croissance[c(2,9)]
names(taux_croissanceF)[names(taux_croissanceF) == "Country.Code"] <- "ISO3"
names(taux_croissanceF)[names(taux_croissanceF) == "X2019"] <- "taux_croissance"

Taux_pop_active <- read.csv("~/Documents/OpenClassRooms/Projet9/Data/Taux_pop_active.csv", encoding="UTF-8", sep=";")
Taux_pop_activeF = Taux_pop_active[c(2,9)]
names(Taux_pop_activeF)[names(Taux_pop_activeF) == "Country.Code"] <- "ISO3"
names(Taux_pop_activeF)[names(Taux_pop_activeF) == "X2019"] <- "Taux_pop_active"

I.3) Nettoyage et création du jeu de donnée final

1) Création d’un DF pour l’importation, l’exportation, la consommation totale et la production de poulet

Création d’une colonne de Consommation totale , égale à la somme de la colonne Food (nourriture), Processing (traitement) et Tourist consumption (consommation par les touristes)

#Format long to format wide by Element (et retrait du Element.code) :
FB2019 = spread(data=FB2019[-c(3,5,8)], key=Element, Value)

#Somme de la colonne Food, Processing et Tourist Consumption
FB2019=FB2019 %>% 
  rowwise() %>% 
  mutate(Consommation_Tot = sum(Food, Processing,`Tourist consumption` , na.rm = TRUE))

#Somme de la colonne Food, Processing et Tourist Consumption
FB2019 = FB2019[, -c(6,8,10)]

#Changement en format Long. 
FB2019 = gather(FB2019, Element, Value,`Export Quantity`:Consommation_Tot, factor_key=TRUE)

#Changement en format Wide, avec une colonne par Item. 
FB2019 <- spread(FB2019, key=Item, Value)

#Changement de nom des colonnes pour plus de clarté
names(FB2019)[5:9] <- c("Cereales","Poisson","Viande","Poulet",'Legumes')

## Création d'un DF pour l'importation, l'exportation, la consommation et la production
Importation = subset(FB2019, Element == 'Import Quantity')
Exportation = subset(FB2019, Element == 'Export Quantity')
Consommation = subset(FB2019, Element == 'Consommation_Tot')
Production = subset(FB2019, Element == 'Production')


## Changement du nom de chaque item pour ajouter son Element. 

for (i in c(5:9)) {
  names(Importation)[i] <- paste(names(Importation)[i],'Importation',sep="_")
  names(Exportation)[i] <- paste(names(Exportation)[i],'Exportation',sep="_")
  names(Consommation)[i] <- paste(names(Consommation)[i],'Consommation',sep="_")
  names(Production)[i] <- paste(names(Production)[i],'Production',sep="_")
}

2) Jointure des données

Jointure des données des indicateurs PESTEL et des indicateurs relatif au commerce et la consommation de poulet. Les jointures sont entières (‘Full Join’, où l’on garde toutes les lignes des jeux de données)

Data_m=merge(PIB2019[,c(1,2,9)],PoliticalStability2019[,c(1,9)], by='ISO3',all = T) 
Data_m=merge(Data_m,Population2019[,c(1,9)], by='ISO3',all = T) 
Data_m=merge(Data_m,Acces_electriciteF, by='ISO3',all = T)
Data_m=merge(Data_m,Emission_CO2F, by='ISO3',all = T)
Data_m=merge(Data_m,Indice_juridiqueF, by='ISO3',all = T)
Data_m=merge(Data_m,Indice_perf_logistiqueF, by='ISO3',all = T)
Data_m=merge(Data_m,Investissement_etrangerF, by='ISO3',all = T)
Data_m=merge(Data_m,taux_croissanceF, by='ISO3',all = T) 
Data_m=merge(Data_m,Taux_pop_activeF, by='ISO3',all = T) 
Data_m=merge(Data_m,Consommation[,c(1,5:9)], by='ISO3',all = T) 
Data_m=merge(Data_m,Exportation[,c(1,5:9)], by='ISO3',all = T) 
Data_m=merge(Data_m,Importation[,c(1,5:9)], by='ISO3',all = T) 
Data_m=merge(Data_m,Production[,c(1,5:9)], by='ISO3',all = T)

## Changement des noms du PIB, Stabilité et Population, modifiés après jointures.  

names(Data_m)[3:5] <- c("PIB","Stabilite","Population")

3) Gestion des valeurs manquantes

## Retrait des colonnes concernant la viande, les céréales, les légumes et le poisson. 

Data_clustering <- Data_m[,c(1,3:12,16,21,26,31)]


## Affichage du nombre de ligne sans valeur

colSums(is.na(Data_clustering))

##                    ISO3                     PIB               Stabilite 
##                       0                      60                      73 
##              Population       Acces_electricite            Emission_CO2 
##                      91                      10                      34 
##        Indice_juridique         Perf_logistique Investissement_etranger 
##                      35                      59                      37 
##         taux_croissance         Taux_pop_active     Poulet_Consommation 
##                      23                      35                      91 
##      Poulet_Exportation      Poulet_Importation       Poulet_Production 
##                     134                      92                      93

# L'information capitale de cette analyse concerne le commerce de poulet. Il est nécessaire de retirer les pays dont les données ne sont pas disponible. 

Data_clustering = Data_clustering %>% drop_na(Poulet_Consommation,Poulet_Importation,Poulet_Production)

colSums(is.na(Data_clustering))

##                    ISO3                     PIB               Stabilite 
##                       0                       1                       3 
##              Population       Acces_electricite            Emission_CO2 
##                       0                       1                       5 
##        Indice_juridique         Perf_logistique Investissement_etranger 
##                       7                      21                       6 
##         taux_croissance         Taux_pop_active     Poulet_Consommation 
##                       5                       7                       0 
##      Poulet_Exportation      Poulet_Importation       Poulet_Production 
##                      40                       0                       0

#Les données du PIB pour la Chine et de la stabilité politique de Taiwan sont trouvable sur d'autres sources de donnée. 

dplyr::filter(Data_clustering,is.na(PIB)) 

Data_clustering$PIB[Data_clustering$ISO3=='TWN']<-25903 # source : countryeconomy

dplyr::filter(Data_clustering,is.na(Stabilite)) 

Data_clustering$Stabilite[Data_clustering$ISO3=='CHN']<- -0.26 # source : Banque Mondiale

#retrait de l'exportation de poulet qui compte 40 pays sans valeurs.
Data_clustering = Data_clustering[-13]



colSums(is.na(Data_clustering))

##                    ISO3                     PIB               Stabilite 
##                       0                       0                       2 
##              Population       Acces_electricite            Emission_CO2 
##                       0                       1                       5 
##        Indice_juridique         Perf_logistique Investissement_etranger 
##                       7                      21                       6 
##         taux_croissance         Taux_pop_active     Poulet_Consommation 
##                       5                       7                       0 
##      Poulet_Importation       Poulet_Production 
##                       0                       0

#retrait pays qui ont au moins une valeur non renseignée. 

Data_clustering = Data_clustering %>% drop_na()

II) Analyse de la relation entre les variables : L’Analyse en Composante Principale

II.1) Mesure de la corrélation entre les variables

#Attribution du code ISO3 en tant qu'ID 
rownames(Data_clustering) = Data_clustering[,1]

#Calcul de la corrélation entre les variables
res <- cor(Data_clustering[-1])
res=round(res, 3) 


#représentation graphique :

corrplot(res, type = "upper", order = "hclust", 
         tl.col = "black", tl.srt = 45)

#La consommation et la production de Poulet sont extrêmement corrélées (r=0.99). Pour simplifier nous ne garderons que la consommation de poulet. 
Data_Clus.def <- Data_clustering[-c(1,14)]

Voici la structure du jeu de donnée final :

head(Data_Clus.def,5)

II.2) Analyse en composante principale

Le but de l’ACP est de synthétiser, réduire le nombre de variables et d’observer les groupes de corrélation entre elles. Cette méthode permet de révéler des groupes d’observations (ici des pays), qui ne pourraient être vus en analysant les variables séparemment.

# Cette commande réalise l'ACP. L'argument scale.unit permet de centrer et réduire les données. 
res.pca <- PCA(Data_Clus.def, scale.unit = TRUE, graph = F, ncp=4)

#Graphe montrant le pourcentage de variance expliqué par chaque axe (dimension)
fviz_eig(res.pca, addlabels = TRUE,main = "", ylab = "Pourcentage d'inertie (variance expliquée)")

Afin de déterminer le nombre de dimensions (axes) à analyser par la suite, nous pouvons garder les dimensions qui expliquent plus de 100/p de la variance, où p=nombre de variables. Ici, nous garderons les 4 premières dimensions qui représentent 67% de la variance totale.

Visualisation des variables sur le cercle des Corrélations

p12=fviz_pca_var(res.pca,
                 ggtheme = theme_minimal(),labelsize = 3, title = "ACP - Cercle des corrélations")
p13=fviz_pca_var(res.pca,
                 ggtheme = theme_minimal(),
                 axes = c(1, 3),labelsize = 3,title = "")
p14=fviz_pca_var(res.pca,
                 ggtheme = theme_minimal(),
                 axes = c(1, 4),labelsize = 3,title = "")

plot_grid(p12,p13 ,p14, labels = "AUTO", nrow=1)

Heatmap de la corrélation, représentation, et contribution des variables sur les axes de l’ACP

#Préparation des couleurs pour les heatmaps
col <- viridis(n = 100,begin = 0,end = 1)

#Représentation des Heatmaps
heatmap.2(as.matrix(res.pca$var$cor), trace="none", main="Corrélation entre les variables \net les composantes principales",
    xlab="Composantes principales",
    ylab="Variables",
    margins = c(6, 16),
    col = col,
    Rowv = F,
    Colv=F,
    dendrogram = 'none',
    cexCol = 1.4)

heatmap.2(as.matrix(res.pca$var$cos2), trace="none", main="Représentation des variables \n dans les composantes principales",
    xlab="Composantes principales",
    ylab="Variables",
    margins = c(6, 16),
    col = col,
    Rowv = F,
    Colv=F,
    dendrogram = 'none',
    cexCol = 1.4)

heatmap.2(as.matrix(res.pca$var$contrib), trace="none", main="Contribution des variables \nsur les composantes principales",
    xlab="Composantes principales",
    ylab="Variables",
    margins = c(6, 16),
    col = col, 
    Rowv = F,
    Colv=F,
    dendrogram = 'none',
    cexCol = 1.4)

Les heatmaps ci-dessus représentent la contribution, la représentation et la corrélation entre les différentes variables et les dimensions de l’ACP. L’importation de poulet, la variable d’intérêt principale, semble bien corrélée et représentée sur l’axe 1 et 2 de l’ACP.

Visualisation des individus sur les principaux axes

p12=fviz_pca_ind(res.pca,
                 ggtheme = theme_minimal(),title = "ACP - Individus")
p13=fviz_pca_ind(res.pca,
                 ggtheme = theme_minimal(),
                 axes = c(1, 3),title = "")
p14=fviz_pca_ind(res.pca,
                 ggtheme = theme_minimal(),
                 axes = c(1, 4),title = "")

plot_grid(p12,p13 ,p14, labels = "AUTO", nrow=1)

#sauvegarde de la valeur de la variance expliqué par chaque axe.
variance.percent <- as.data.frame(round(res.pca$eig, 2))

#sauvegarde de la position des individus dans les coordonnées de l'ACP. 
ind.coord <- as.data.frame(res.pca$ind$coord)

Le graphe ci-dessus, présente la position de chaque pays sur les différents axes de l’ACP. Afin de regrouper les pays relativement similaires entre eux, nous utiliserons deux méthode : la méthode de la Classification Ascendante Hierarchique (ACH) et la méthode des kmeans.

III) Analyse des groupes d’individus

III.1) Regroupement sur les axes de l’ACP.

III.1.a) Classification Ascendante Hiérarchique (CAH)

La première étape est de visualiser le nombre de cluster optimal selon 3 méthodes : gap-stat, méthode des silhouettes et du coude.

p1=fviz_nbclust(ind.coord,hcut, method="gap_stat", verbose = F)+
  geom_vline(xintercept = 6 , linetype = 2, color = 'red')+
  labs(title= "Nombre de clusters optimal : gap statistic") + 
  xlab("Nombre de clusters k")

p2=fviz_nbclust(ind.coord,hcut, method = "silhouette")+
  geom_vline(xintercept = 6 , linetype = 2, color = 'red')+
  labs(title= "Nombre de clusters optimal : méthode des silouhettes") + 
  xlab("Nombre de clusters k")

p3=fviz_nbclust(ind.coord,hcut, method="wss")+
  geom_vline(xintercept = 6 , linetype = 2, color = 'red')+
  labs(title= "Nombre de clusters optimal : méthode du coude") + 
  xlab("Nombre de clusters k")


plot_grid(p1,p2,p3, labels = "AUTO", nrow=1)

Selon ces différents modèles, 6 groupes semblent cohérent.

#Fonction permettant de générer une CAH sur les résultats de l'ACP
res.hcpc <- HCPC(res.pca, graph = F, nb.clust = 7 ,consol = F )

#sauvegarde de la valeur de chaque groupe attribué aux individus
ind.coord$CAH_clust <- res.hcpc$data.clust$clust

#représentation graphique : Dendrogramme puis visualisation des individus selon les 3 dim majoritaires de l'ACP. 
p.dendo=fviz_dend(res.hcpc, cex = 0.7, k =7,palette = "jco", color_labels_by_k = T,rect = TRUE, rect_fill = TRUE, 
          rect_border = "jco", main = "Dendrogramme" ,ylab = "Distance entre clusters (méthode de Ward)",xlab = "Pays")

p.CAH.12=ggscatter(
  ind.coord, x = "Dim.1", y = "Dim.2", 
  color = "CAH_clust", palette =c('#0073C2FF','#EFC000FF',"#868686FF",'#7AA6DCFF','#003C67FF','#8F7700FF','#CD534CFF'),
  ellipse = TRUE, ellipse.type = "convex", ellipse.alpha = 0.15,
  size = 1.5,  legend = "right", ggtheme = theme_minimal(),
  xlab = paste0("Dim 1 (", variance.percent[1,2], "% )" ),
  ylab = paste0("Dim 2 (", variance.percent[2,2], "% )" ),
  ) +
  geom_text_repel(aes(label = rownames(ind.coord),color=CAH_clust,fontface="bold"),size = 3.5)

p.CAH.13=ggscatter(
  ind.coord, x = "Dim.1", y = "Dim.3", 
  color = "CAH_clust", palette =c('#0073C2FF','#EFC000FF',"#868686FF",'#7AA6DCFF','#003C67FF','#8F7700FF','#CD534CFF'),
  ellipse = TRUE, ellipse.type = "convex", ellipse.alpha = 0.15,
  size = 1.5,  legend = "right", ggtheme = theme_minimal(),
  xlab = paste0("Dim 1 (", variance.percent[1,2], "% )" ),
  ylab = paste0("Dim 3 (", variance.percent[3,2], "% )" ),
  ) +
  geom_text_repel(aes(label = rownames(ind.coord),color=CAH_clust,fontface="bold"),size = 3.5)

p.CAH.14=ggscatter(
  ind.coord, x = "Dim.2", y = "Dim.3", 
  color = "CAH_clust", palette =c('#0073C2FF','#EFC000FF',"#868686FF",'#7AA6DCFF','#003C67FF','#8F7700FF','#CD534CFF'),
  ellipse = TRUE, ellipse.type = "convex", ellipse.alpha = 0.15,
  size = 1.5,  legend = "right", ggtheme = theme_minimal(),
  xlab = paste0("Dim 1 (", variance.percent[1,2], "% )" ),
  ylab = paste0("Dim 4 (", variance.percent[4,2], "% )" ),
  ) +
  geom_text_repel(aes(label = rownames(ind.coord),color=CAH_clust,fontface="bold"),size = 3.5)


bottom_row <- plot_grid(plot_grid(p.CAH.12, p.CAH.13, p.CAH.14, labels = c('B', 'C','D'),ncol = 3))
                        
plot_grid(p.dendo, bottom_row, labels = c('A', ''), ncol = 1, rel_heights = c(.9, 1.2))

heatmaply(t(as.matrix(ind.coord[1:4])), hclust_method = "ward.D2",fontsize_col = 6, k_col = 7,Rowv = F)

Le groupe n°6 dispose des plus hautes valeurs au niveau de la dimension 1, où l’importation de poulet est bien représentée.

# Ajout des groupes aux données d'origine
Data_Clus.def$ACP_CAH <- as.factor(res.hcpc$data.clust$clust)
subset(rownames(Data_Clus.def), Data_Clus.def$ACP_CAH == 6 )

##  [1] "ARE" "AUS" "AUT" "BEL" "CAN" "CHE" "CZE" "DEU" "DNK" "ESP" "EST" "FIN"
## [13] "FRA" "GBR" "IRL" "ISL" "JPN" "KOR" "LUX" "NLD" "NOR" "NZL" "SWE"

III.1.b) Regroupement par la méthode des kmeans

#Visualisation du nombre de cluster optimal selon 3 méthodes
p1=fviz_nbclust(ind.coord[1:4],kmeans, method="gap_stat", verbose = F)+
  geom_vline(xintercept = 6 , linetype = 2, color = 'red')+
  labs(title= "Nombre de clusters optimal : gap statistic") + 
  xlab("Nombre de clusters k")

p2=fviz_nbclust(ind.coord[1:4],kmeans, method = "silhouette")+
  geom_vline(xintercept = 6 , linetype = 2, color = 'red')+
  labs(title= "Nombre de clusters optimal : méthode des silouhettes") + 
  xlab("Nombre de clusters k")

p3=fviz_nbclust(ind.coord[1:4],kmeans, method="wss")+
  geom_vline(xintercept = 6 , linetype = 2, color = 'red')+
  labs(title= "Nombre de clusters optimal : méthode du coude") + 
  xlab("Nombre de clusters k")


plot_grid(p1,p2,p3, labels = "AUTO", nrow=1)

#Selon la méthode du coude (wss) et minoritairement de la méthode des silhouettes et gap stat, 6 groupes semblent cohérent. 

km.out <- kmeans(ind.coord[1:4], centers = 6)


#sauvegarde de la valeur de chaque groupe attribué aux individus
ind.coord$kmeans_clust <- factor(km.out$cluster)

#représentation graphique 
km.out$centers

##        Dim.1       Dim.2      Dim.3      Dim.4
## 1 -1.3863926  0.43189292 -0.8872394 -0.9798607
## 2 -0.2185248 -0.35567577  1.1805305 -0.3228578
## 3  1.2460585 -0.09200155 -0.5166095 -0.2435474
## 4  3.2851417 -0.71826315 -0.5263452  0.4537405
## 5 -2.0351942 -0.09428507 -0.1945350  1.3780941
## 6  3.2911422  7.24485183  1.7533678  0.6498309

heatmap.2(t(as.matrix(km.out$centers)), trace="none", main="Heatmap de la position des centroids \n des clusters sur les différentes dimensions de l'ACP",
    xlab="Clusters (kmeans)",
    ylab="Variables",
    margins = c(8, 8),
    col = col,
    Rowv = F,
    Colv=F,
    dendrogram = 'none')

Sur les dimensions 1 et 2, Le groupe numéro 4 et 6 ont les plus forte valeur sur les dimensions 1 et 2. Le groupe 6 cependant est largement biaisé par la population de ces pays très importante (Chine, USA, IND).

pkmeans12=ggscatter(
  ind.coord, x = "Dim.1", y = "Dim.2", 
  color = "kmeans_clust", palette = "jco",
  ellipse = TRUE, ellipse.type = "convex", ellipse.alpha = 0.15,
  size = 1.5,  legend = "right", ggtheme = theme_minimal(),
  xlab = paste0("Dim 1 (", variance.percent[1,2], "% )" ),
  ylab = paste0("Dim 2 (", variance.percent[2,2], "% )" ),
  ) +
  stat_mean(aes(color = kmeans_clust), size = 4)+
  geom_text_repel(aes(label = rownames(ind.coord),color=kmeans_clust,fontface="bold"),size = 3.5)
 
pkmeans13=ggscatter(
  ind.coord, x = "Dim.1", y = "Dim.3", 
  color = "kmeans_clust", palette = "jco",
  ellipse = TRUE, ellipse.type = "convex", ellipse.alpha = 0.15,
  size = 1.5,  legend = "right", ggtheme = theme_minimal(),
  xlab = paste0("Dim 1 (", variance.percent[1,2], "% )" ),
  ylab = paste0("Dim 3 (", variance.percent[3,2], "% )" ),
  ) +
  stat_mean(aes(color = kmeans_clust), size = 4)+
  geom_text_repel(aes(label = rownames(ind.coord),color=kmeans_clust,fontface="bold"),size = 3.5)

pkmeans14=ggscatter(
  ind.coord, x = "Dim.1", y = "Dim.4", 
  color = "kmeans_clust", palette = "jco",
  ellipse = TRUE, ellipse.type = "convex", ellipse.alpha = 0.15,
  size = 1.5,  legend = "right", ggtheme = theme_minimal(),
  xlab = paste0("Dim 1 (", variance.percent[1,2], "% )" ),
  ylab = paste0("Dim 4 (", variance.percent[4,2], "% )" ),
  ) +
  stat_mean(aes(color = kmeans_clust), size = 4)+
  geom_text_repel(aes(label = rownames(ind.coord),color=kmeans_clust,fontface="bold"),size = 3.5)

plot_grid(pkmeans12,pkmeans13,pkmeans14, labels = "AUTO", nrow=1)

# Ajout des groupes aux données d'origine

Data_Clus.def$ACP_kmeans <- as.factor(km.out$cluster)

subset(rownames(Data_Clus.def), Data_Clus.def$ACP_kmeans == 4)

##  [1] "ARE" "AUS" "AUT" "BEL" "CAN" "CHE" "CZE" "DEU" "DNK" "FIN" "FRA" "GBR"
## [13] "IRL" "ISL" "JPN" "LUX" "NLD" "NOR" "NZL" "SWE"

III.2) Regroupement sur les variables de départ

III.2.a) Classification Ascendante Hiérarchique (CAH)

## première étape : centrer et réduire les données
Data_clustering.scale=scale(Data_Clus.def[1:12])

## Deuxième étape : calculer la distance (euclidienne) entre les points
Data_clustering.scal.dist=dist(Data_clustering.scale)

## Classification hiérarchique : 
hc.out <- hclust(Data_clustering.scal.dist, method = "ward.D2")

## détermination du nombre de groupe nécessaire: 
p1=fviz_nbclust(Data_clustering.scale,hcut, method="gap_stat", verbose = F)+
  geom_vline(xintercept = 7 , linetype = 2, color = 'red')+
  labs(title= "Nombre de clusters optimal : gap statistic") + 
  xlab("Nombre de clusters k")

p2=fviz_nbclust(Data_clustering.scale,hcut, method = "silhouette")+
  geom_vline(xintercept = 7 , linetype = 2, color = 'red')+
  labs(title= "Nombre de clusters optimal : méthode des silouhettes") + 
  xlab("Nombre de clusters k")

p3=fviz_nbclust(Data_clustering.scale,hcut, method="wss")+
  geom_vline(xintercept = 7 , linetype = 2, color = 'red')+
  labs(title= "Nombre de clusters optimal : méthode du coude") + 
  xlab("Nombre de clusters k")

plot_grid(p1, p2, p3, labels = "AUTO",ncol = 3)

## Selon la méthode du coude et des silhouette : 7 groupes semblent cohérents

pvizdend=fviz_dend(hc.out, cex = 0.6, k =7,palette = "jco",rect = TRUE, color_labels_by_k = F, rect_border ="jco", rect_fill = TRUE, lower_rect = -3)

## Cluster
hc.out.cluster7 <- cutree(hc.out,k = 7)


p12=fviz_cluster(list(data=Data_clustering.scale, cluster = hc.out.cluster7), repel = TRUE, 
                 show.clust.cent = FALSE, 
                 axes = c(1, 2),
                 palette = c("#EFC000FF",'#0073C2FF','#8F7700FF',"#7AA6DCFF",'#CD534CFF',"#3B3B3BFF",'#003C67FF'), 
                 ggtheme = theme_minimal(),
                 main="Clustering sur les variables : CAH")

p13=fviz_cluster(list(data=Data_clustering.scale, cluster = hc.out.cluster7), repel = TRUE, 
                 show.clust.cent = FALSE, 
                 axes = c(1, 3),
                 palette =c("#EFC000FF",'#0073C2FF','#8F7700FF',"#7AA6DCFF",'#CD534CFF',"#3B3B3BFF",'#003C67FF'),
                 ggtheme = theme_minimal(),
                 main="")

p14=fviz_cluster(list(data=Data_clustering.scale, cluster = hc.out.cluster7), repel = TRUE, 
                 show.clust.cent = FALSE, 
                 axes = c(1, 4),
                 palette = c("#EFC000FF",'#0073C2FF','#8F7700FF',"#7AA6DCFF",'#CD534CFF',"#3B3B3BFF",'#003C67FF'), 
                 ggtheme = theme_minimal(),
                 main="")

bottom_row <- plot_grid(plot_grid(p12, p13, p14, labels = c('B', 'C','D'),ncol = 3))
                        
plot_grid(pvizdend, bottom_row, labels = c('A', ''), ncol = 1, rel_heights = c(1, 1.2))

heatmaply(t(Data_clustering.scale), hclust_method = "ward.D2",fontsize_col = 6, k_col = 7,Rowv = F)

Le groupe 4 dispose de la plus forte importation de Poulet et des indicateurs socio-économiques relativement important.

Data_Clus.def$CAH <- as.factor(hc.out.cluster7)

subset(rownames(Data_Clus.def), Data_Clus.def$CAH == 4)

##  [1] "ARE" "BEL" "DEU" "FRA" "GBR" "JPN" "MEX" "NLD" "SAU" "ZAF"

III.2.b) Regroupement par la méthode des kmeans

## Nombre de clusters
p1=fviz_nbclust(Data_clustering.scale,kmeans, method="gap_stat", verbose = F)+
  geom_vline(xintercept = 5 , linetype = 2, color = 'red')+
  labs(title= "Nombre de clusters optimal : gap statistic") + 
  xlab("Nombre de clusters k")

p2=fviz_nbclust(Data_clustering.scale,kmeans, method = "silhouette")+
  geom_vline(xintercept = 5 , linetype = 2, color = 'red')+
  labs(title= "Nombre de clusters optimal : méthode des silouhettes") + 
  xlab("Nombre de clusters k")

p3=fviz_nbclust(Data_clustering.scale,kmeans, method="wss")+
  geom_vline(xintercept = 5 , linetype = 2, color = 'red')+
  labs(title= "Nombre de clusters optimal : méthode du coude") + 
  xlab("Nombre de clusters k")
plot_grid(p1, p2, p3, labels = "AUTO", nrow = 1)

## Selon la méthode des silhouettes et du gap statistique: 5 groupes semblent cohérents

km.out.2 <- kmeans(Data_clustering.scale, centers = 5, nstart = 100)

#sauvegarde de la valeur de chaque groupe attribué aux individus

km.cluster2 <- km.out.2$cluster

#représentation graphique 

col2 <- viridis(n = 100)


heatmap.2(t(as.matrix(km.out.2$centers)), trace="none", main="Heatmap de la position des centroids \n sur les différentes variables",
    ylab="Variables",
    xlab="Clusters (kmeans)",
    margins = c(6, 15),
    col = col2,
    Rowv = F,
    Colv=F,
    dendrogram = 'none')

Le groupe numéro 1 et 3 semble avoir une forte importation de Poulet. Pour Le groupe n°3, l’importation de poulet (et la consommation) est probablement liée à une forte population.

pkmeans12=fviz_cluster(list(data=Data_clustering.scale, cluster = km.cluster2), repel = TRUE, 
                       show.clust.cent = T, 
                       axes = c(1, 2),
                       palette = "jco", 
                       ggtheme = theme_minimal(), 
                       main = "Factor map")

pkmeans13=fviz_cluster(list(data=Data_clustering.scale, cluster = km.cluster2), repel = TRUE, 
                       show.clust.cent = T, 
                       axes = c(1, 3),
                       palette = "jco", 
                       ggtheme = theme_minimal(), 
                       main = "Factor map")
pkmeans14=fviz_cluster(list(data=Data_clustering.scale, cluster = km.cluster2), repel = TRUE, 
                       show.clust.cent = T, 
                       axes = c(1, 4),
                       palette = "jco", 
                       ggtheme = theme_minimal(), 
                       main = "Factor map")


plot_grid(pkmeans12, pkmeans13, pkmeans14, labels = "AUTO", nrow=1)

Data_Clus.def$kmeans <- as.factor(km.cluster2)

subset(rownames(Data_Clus.def), Data_Clus.def$kmeans == 1)

##  [1] "ARE" "AUS" "AUT" "BEL" "CAN" "CHE" "CZE" "DEU" "DNK" "ESP" "EST" "FIN"
## [13] "FRA" "GBR" "IRL" "ISL" "ISR" "ITA" "JPN" "KOR" "KWT" "LUX" "MEX" "NLD"
## [25] "NOR" "NZL" "OMN" "PRT" "SAU" "SVN" "SWE" "ZAF"

Bilan

# Remise du code ISO3 en variable et rajout du nom des pays
Data_Clus.def.ISO<-Data_Clus.def

Data_Clus.def.ISO$ISO3 <- rownames(Data_Clus.def.ISO)
rownames(Data_Clus.def.ISO) <- NULL

Data_Clus.def.ISO <- left_join(Data_Clus.def.ISO,Data_m[c(1,2)],by='ISO3')

Reduce(intersect, list(subset(Data_Clus.def.ISO[c(14,18)], ACP_kmeans == 4)$Area,
                       subset(Data_Clus.def.ISO[c(13,18)], ACP_CAH == 6)$Area,
                       subset(Data_Clus.def.ISO[c(16,18)], kmeans == 1)$Area,
                       subset(Data_Clus.def.ISO[c(15,18)], CAH == 4)$Area))

## [1] "United Arab Emirates"                                
## [2] "Belgium"                                             
## [3] "Germany"                                             
## [4] "France"                                              
## [5] "United Kingdom of Great Britain and Northern Ireland"
## [6] "Japan"                                               
## [7] "Netherlands"

Data_Clus.def.ISO$Groupe_Focus <- ifelse(Data_Clus.def.ISO$ISO3 %in% c("BEL","FRA","JPN","ARE","DEU","GBR","NLD"),'1', '0')

Globalement, on distingue sept pays communs dans les quatre méthodes de regroupement. Ces pays ont une importante importation de Poulet et des bonnes conditions économiques et sociales. La France fait également partie de cette liste de pays cible, où l’entreprise La Poule qui Chante s’est bien développée. Nous pouvons donc valider cette liste de pays cible.

Exportation en Europe : Belgique, Pays-Bas, Allemagne, Royaume Uni.

Exportation hors de l’Europe : Japon et Emirat Arabe Unis.

“FRA”,“BEL”,“JPN”,“ARE”,“DEU”,“GBR”,“NLD”

Résumé des indicateurs des pays cibles

Data_selected = subset(Data_Clus.def.ISO, ISO3 %in% c("BEL","FRA","JPN","ARE","DEU","GBR","NLD") )
Data_selected$Area[Data_selected$ISO3=='GBR']<-'UK'
Data_selected2 = gather(Data_selected, Element, Value,`PIB`:Poulet_Importation, factor_key=TRUE)

ggplot(data= Data_selected2, aes(x=reorder(Area,Value), y=Value))+
  geom_bar(stat="identity",width = .8)+
 # facet_grid(Element~.,scales = "free")+
  facet_wrap(.~Element,scales = "free",ncol=4)+
  labs( y="Valeurs des indicateurs", x="Pays" )+
  theme(axis.text.x = element_text(angle = 45, hjust=1))+ coord_flip()

Carte des pays cibles

my_spdf <- readOGR("~/Documents/OpenClassRooms/Projet9/Data/TM_WORLD_BORDERS_SIMPL-0.3/TM_WORLD_BORDERS_SIMPL-0.3.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/dubois/Documents/Openclassrooms/Projet9/Data/TM_WORLD_BORDERS_SIMPL-0.3/TM_WORLD_BORDERS_SIMPL-0.3.shp", layer: "TM_WORLD_BORDERS_SIMPL-0.3"
## with 246 features
## It has 11 fields
## Integer64 fields read as strings:  POP2005

DATA.MAP <- merge(my_spdf,Data_Clus.def.ISO, by="ISO3")



# Create a color palette for the map:

pal <- colorFactor(
  palette = c("#DDCC77","#D55E00"),
  domain = DATA.MAP@data$Groupe_Focus
)

# Prepare the text for tooltips:

mytext <- paste(
   "Pays: ", DATA.MAP@data$NAME,"<br/>", 
   "Stabilité politique: ", DATA.MAP@data$Stabilite,"<br/>",
   "Population: ", DATA.MAP@data$Population,"<br/>",
   "Performance logistique: ", round(DATA.MAP@data$Perf_logistique,2),"<br/>",
   "Investissements étrangers directs (% du PIB): ",round(DATA.MAP@data$Investissement_etranger,2),"<br/>",
   "Consommation de Poulet (t): ", DATA.MAP@data$Poulet_Consommation,"<br/>",
   "Importation de Poulet (t): ", DATA.MAP@data$Poulet_Importation,"<br/>",
   "ACP_CAH: ", DATA.MAP@data$ACP_CAH,"<br/>",
   "ACP_kmeans: ", DATA.MAP@data$ACP_kmeans,"<br/>",
   "CAH: ", DATA.MAP@data$CAH,"<br/>",
   "Kmeans: ", DATA.MAP@data$kmeans,"<br/>",
  sep="") %>%
  lapply(htmltools::HTML)


mydf <- data.frame(Observation = c("FRA","BEL","JPN","ARE","DEU","GBR","NLD"),
                   InitialLat = subset(DATA.MAP, ISO3 %in%  c("FRA"))$LAT,
                   InitialLong = subset(DATA.MAP, ISO3 %in%  c("FRA"))$LON,
                   NewLat = subset(DATA.MAP, ISO3 %in%  c("FRA","BEL","JPN","ARE","DEU","GBR","NLD"))$LAT,
                   NewLong =subset(DATA.MAP, ISO3 %in%   c("FRA","BEL","JPN","ARE","DEU","GBR","NLD"))$LON,
                   stringsAsFactors = FALSE)

subset(DATA.MAP, ISO3 %in%  c("BEL","FRA","JPN","ARE","DEU","GBR","NLD"))$LON

## [1]   2.550   9.851 139.068   4.664   5.389  -1.600  54.163

subset(DATA.MAP, ISO3 %in%  c("BEL"))$LAT

## [1] 50.643

mydf2 <- data.frame(group =  c("FRA","BEL","JPN","ARE","DEU","GBR","NLD"),
                    lat = subset(DATA.MAP, ISO3 %in%   c("FRA","BEL","JPN","ARE","DEU","GBR","NLD"))$LAT,
                    long = subset(DATA.MAP, ISO3 %in%  c("BEL","FRA","JPN","ARE","DEU","GBR","NLD"))$LON)

mydf2 <- data.frame(group = c("FRA","BEL","JPN","ARE","DEU","GBR","NLD"),
                    lat = c(mydf$InitialLat, mydf$NewLat),
                    long = c(mydf$InitialLong, mydf$NewLong))

# Final Map
m <- leaflet(DATA.MAP) %>% 
  addTiles()  %>% 
  setView( lat=10, lng=5 , zoom=2) %>%
  addPolygons( 
    fillColor = ~pal(DATA.MAP@data$Groupe_Focus), 
    stroke=TRUE, 
    fillOpacity = 0.9, 
    color="black", 
    weight=0.3,
    label = mytext,
    labelOptions = labelOptions( 
      style = list("font-weight" = "normal", padding = "3px 8px"), 
      textsize = "13px", 
      direction = "auto"
    )
  ) %>%
  addLegend( pal=pal, values=~DATA.MAP@data$Groupe_Focus, opacity=0.9, title = "Pays cibles", position = "bottomleft" )%>%
addPolylines(data = mydf2[c(1,2,8,9),], lng = ~long, lat = ~lat, group = ~group,weight = 2)%>%
addPolylines(data = mydf2[c(1,3,8,10),], lng = ~long, lat = ~lat, group = ~group,weight = 2)%>%
addPolylines(data = mydf2[c(1,4,8,11),], lng = ~long, lat = ~lat, group = ~group,weight = 2)%>%
addPolylines(data = mydf2[c(1,5,8,12),], lng = ~long, lat = ~lat, group = ~group,weight = 2)%>%
addPolylines(data = mydf2[c(1,6,8,13),], lng = ~long, lat = ~lat, group = ~group,weight = 2)%>%
addPolylines(data = mydf2[c(1,7,8,14),], lng = ~long, lat = ~lat, group = ~group,weight = 2)

m